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Abstract 

We present a method of adaptive 
grid refinement for the solution of the 
steady Euler equations for transonic flow. 
Our algorithm automatically decides 
where the coarse grid accuracy is insuffi- 
cient, and creates locally uniform refined 
grids in these regions. Tliis typically 
Occurs at the leading and trailing edges. 
The solution is then integrated to steady 
state using the same integrator (FL052) 
in the interior of cr»ch grid. We examine 
the boundary conditions needed on the 
fine grids, and discuss the importance of 
treating the finc/coarse grid interface con- 
servatively. Numeric^ results arc 
presented. 

1. Introduction 

In computing transonic flow fields 
about complex geometries, it is difficult to 
resolve all features of the solution to the 
same accuracy with a uniform grid. As 
much as possible, the regions where the 
solution needs finer grid resolution are 
finely zoned in the initial (pre-solution) 
grid generation phase. However, it is not 
^ways known in advance where those 
regions arc, or how finely zoned to make 
them. The location of the inaccurate 
regions changes vrith different flow 
parameters, mach number, angle of 
attack, etc. 

“Supported in part by Department of Energy Qxi- 
tract No, DEAC0276ER03077-V. 

+ Supported in part by the CXfice of Naval 
Research under Q"ant N00014-81-K-0379 and by 
NASA Langley Researdi Center under Grant 
NAG* 1-186, 


Algorithms arc commonly found in 
the literature where the user computes a 
solution, le-grids, and re-solvcs [1]. In 
this paper, we present an algorithm for 
autoraatir local grid refinement We 
describe a simple prooc^duTe to discover 
the regions of high error (typically the 
leading and trailing edges and in the 
neighborhood of shodc waves), and to 
re-grid by introducing any number of 
lo(^ rectangular fine grids. This both 
removes the guesswork and obtains com- 
parable solutions at less cost than those 
obtained by uniformly refining the grid 
over the entire flow field. 

A wide variety of approaches to 
adapting the grid for better solution reso- 
lution have bren tried. Rai and Ander- 
son [2] have a method of clustering the 
grid lines in the neighborhood of a ^ock 
by "attracting" the lines into the region. 
Harten and Hyman [3] have an algorithm 
where each grid point can move within a 
base grid ocU which stays fixed. In one 
dimension this method can have the same 
sharp resolution as a shock fitting 
scheme. Recent work by Usab and Mur- 
man [4] proposed grid refinement pro- 
cedures similar to ours, but does not 
incorporate the automatic error estimation 
in our approach. 

In section 2 of this paper, we 
describe the algorithm for local grid 
refinement, that is, the error estimation 
and grid generation. We also describe 
how to integrate the solution on these 
multiple grids to steady state, which we 
do using FL052 [5]. Since the refined 
grids arc locally uniform patches in the 
same coordinate system as the coairte 



» 2 “ 


grid, wc arc able to use an existing 
integration routine with very little modifi- 
cation, Section 3 deals with the boundary 
conditions needed on the fine grid. Our 
strategy is to solve an initial boundary 
value problem on each grid. If a fine grid 
touches the airfoil, or has a farficld boun- 
dary, the seme boundary conditions arc 
applied as would be used for a single grid 
computation. The only new type of 
boundary that arises is the finaooarsc 
grid interface. We discuss the importance 
of treating this interface conservatively, 
even if the interface is in a smooth flow 
region, and describe in some detail the 
procedure which wc implement. Finally, 
section 4 compares the multiple grid 
results to a single grid finely zoned run. 

The same issues that arise in the 
interfaces bctv/ccn fine and coarse grids 
(conservation, the data structures and 
bockceping needed for this information), 
arise in the solution of a problem with 
complex geometry by component grids. 
By the latter wc mean multiple grids in 
diJffercnt coordinate systems. In the 
future wc intend to apply our results in 
that direction. Also, since our algorithm 
keeps grids locally uniform, a simple user 
interface is possible. This allows for 
example, the use of a vectorized integra- 
tor. ihe method does not have the major 
drawbacks of moving grid point methods, 
namely, grid skewness and the ”all points 
to the worst zone" problem, and thus 
seems very suitable for 3D calculations as 
well. In this paper, wc present a sys- 
tematic study to verify that this method 
works. Wc demonstrate that with no loss 
in the convergence rate, wc can capture 
the accuracy of the solution on a grid 
twice as fine by using a coarse, global 
grid, and adaptively refining only those 
regions where the error is high. 

2. Multiple Grid Method of Adaptive 
Refinement 

'iliis section presents the algorithm 
wc use to solve the 2D Euler eqi;iations 
for steady flow about an airfoil. We 
describe the overall algorithm before 
going into detail about the main steps. A 
more detailed discussion of the structure 
of this algorithm is found in [6]. 


The solution procedure (described in 
section 2,3) starts by time- stepping on a 
single global grid. Since the initial condi- 
tions are uniform flow, wc wait until the 
solution has settled down to, say, a resi- 
dual ~ 10"^ before applying the error 
estimator and subsequent adaptive stra- 
tegy. The error estimator (described in 
section 2.1) is then applied at every point 
on the coarse grid. Those grid points 
where the estimate is high arc flagged as 
needing finer grid resolution. The grid 
generation algorithm creates fmc grids in 
the same coordinate system as the coarse 
grid, so that every flagged point is con- 
tained in a fine grid. An important point 
is that the fine grids arc rectangles in the 
computational plane. For example, the 
refinement at the leading edge in figure 
2.1a is the center rectangle in the compu- 
tational domain shown in figure 2.1b. 
Since the grid is periodic in the | direc- 
tion, with the break at the trailing edge, 
the trailing edge refined grids are the left- 
most and rightmost rectangles in figure 
2.1b. 




Figure 2.1 Fine grids at the leading 
and trailing edges. 

The use of rectangles as the basis for 
refinement is a crucial decision. First, it 
allows for a very simple user interface. 
The integrator which is used on the global 
coarse grid can also be used without 




" 3 " 


change on each fine grid. Second, the 
use of rectangles makes the data structuie 
problem tractable, since only four comer 
points arc needed to fix the subgrid loca- 
tion. The storage overhead is thus on a 
per grid basis, rather than on a per grid 
point basil, and is negligible. oSjicr 
methods typically use pointers for each 
refined cell of the coarse grid, or possibly 
each row. Finally, this approach to adap- 
tive grid refinement docs not suffer frorc? 
the two main problems in moving grid 
fjoint methods. These problems arc the 
difficulty in controlling grid skewness, 
and the problem of adequately resolving 
several features of the solution when all 
points rush to the strongest feature [71. It 
is clear that some mcd.vnism of adding 
points in a simple way as well as moving 
them is called for. hi our method, if a 
refined grid is found to need further 
refinement (the error estimate at fmc grid 
points is still too high), another, finer rec- 
tangle is added which will be nested In 
the existing subgrid in the same way the 
subgrid is nested in the coarse grid. 

We emphasize that these grids are 
not patched into one global grid, but are 
kept independently, each with its own 
solution vector. This means that some 
coarse grid solution storage is wasted 
(unless it participates in the solution pro- 
cess itself, as in a multigrid method), 
since we will always use the fine grid 
solution when it exists. The benefits 
seem to greatly outweig j this waste, since 
by preventing fragmentation the solution 
process on each grid can stiJi be vector- 
ized, and the loss in computing some 
extra coarse grid points is offset by the 
gain in efficiency due to regularity and 
the simplicity of the data structures. 

Given tMs grid structure, the solu- 
tion on each grid is initialized by interpo- 
lation from the coarse grid, and the time 
stepping continues. Section 2.3 describes 
the integration strategy for multiple grids, 
and reviews both the finite volume 
discretization scheme and the generalized 
Rungc Kutta time stepping method which 
is used to advance the solution on eaeh 
grid. 


2,t, Error Estimation 

In regions of smooth Hew, the cri- 
terion we use for refining the grid is an 
estimate of the error in the solution on 
the finest existing grid in that region. 
Although there is no theory for equations 
of mixed type, in the pmely cUlptie or 
purely hyperbolic case tourc arc estimates 
for the global error in the solution in 
terms of the local truncation error [8]. 
Accordingly, we will estimate the local 
truncation error in the solution using 
ideas similar to Richardson extrapolation 
or deferred correction [6]. To solve 

/( m)jc + S(u)y = 0 
we compute 

fiWr/ = 0, 

where 1/ is the numerical approximation 
to u, and the difference operators approx- 
imating fg and gy based on a stepsize h 
are in O. The lo^ truncation error is 

Q(h)u = rhP, 

where p - 2 for a second order method. 
The term t contains derivatives of the 
solution u. The goid of refining is to 
determine when t is big, and reduce h, 
so that the mme accuracy is attained over 
the entire flow field. T^ idea is to esti- 
mate the error using Richardson 
extrapolation-type estimates by differenc- 
ing on a grid with mesh spacing 2h using 
every otner point of the computed solu- 
tion U. We compute 

Qi2h)U « (2P-1) T hP. 

In the steady state calculation, the resi- 
dual Q{h)U is driven to zero, but the 
coarsen^ grid residual will not be zero. 
Thus we can use 

QWv ^ T HP 

7P-1 

as an estimate of the error at each point. 

Notice that it is unnecessary to know 
the exact form of the truncation error r 
for this method. Also, this residual calcu- 
lation is identical to the first stage of the 
regular Run^e Kutta integration step but 
an a 2h grid. For the error estimation 
step, the computer code is merely 
changed to read i+2 instead of i+1 (and 
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50 on) whon updating the point. (In 
fact, this can somctini» be done automat- 
ically in Fortran by changing the dimen- 
sion statements when dedaring the arrays 
in an integration subroutine). The oora- 
putational overhead of this estimator is 
thus less than one extra integration step 
for the whole computation. 

In regions where the flow is discon- 
tinuous, this procedure no longer gives a 
valid error estimate. However, it acts as 
a trigger for mesh refinement in the pres- 
ence of a Strang shock. We have found 
in our results however, that the largest 
errors arc at the leading edge and ^n 
trailing edge of the airfoil. Of course, this 
depends on the underlying global grid 
resolution. For strong shocks, a fine grid 
is also created in the region of the shock. 
This general procedure has the advantage 
that it is not ..eccssary to know' the loca- 
tion of the shock before the start of the 
computation. 

2.2. Grid Generation 

The output from the error estimation 
routine is a list of (coarse) grid points 
with high error cstinnitcs, indicating that 
a refined zone is needed in that region. 
The grid generation routine separates the 
points into appropriate groups so that a 
(logically) rectangular grid can be placed 
around each group. This proceeds in two 
phases. First, the points that are in dif- 
ferent parts of the domain (such as lead- 
ing edge and trailing edge) arc separated 
into different groups. Around each 
group, a rectangle is formed which is 
large enough to include all points in that 
group. This new grid is then slightly 
enlarged (by one or two coarse grid 
points), to ensure that the fine grid boun- 
dary, where special interpolation formu- 
las will be used for the boundary dif- 
ferencing, is in a region with a small error 
estimate. Figure 2.2 illustrates this pro- 
cedure schematically. 

The only possible exception to this 
procedure is ^own in Figure 2.3. It may 
happen that two different grids are 
created around one cluster of points, in 
order to minimize the size of the 
(unnecessarily) refined region. Details of 
this exception^ ease can te found in [9]. 
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Figure 2.2 Fme grid generation 

around flagged points. 

We make one last remark about the 
rectangular grids. It may happen that tte 
zone needing lefiiKsraent is oblique to the 
underlying grid, for example, foir an 
oblique shock. It could be advantageous 
to be able to align the fmc grid so that 
the coordinates were approximately nor- 
mal and tangent to the di.scontinuity. A 
rotated difference i;cbcmc has been used 
by Jameson [10] for the potential equa- 
tion. Recent results by Davis [11] for the 
Euler equations show much better perfor- 
mance for first order upwind schemes if 
they arc rotated to align with the shock. 
The grid generation procedure has the 
capability to produce grids with this align- 
ment property. However, interpolation 
procedures have not yet been developed 
which treat the fine/coarse grid interface 
conservatively. Work is still in progress 
on this point. 

2.3. Integration Procedure 

It is very easy to solve the equations 
on the grid stn 4 c:turc described above. We 
take one step on all grids using the 
Runge Kutta finite-volume integrator 
^scribed below. Since we ate interested 
in the stciady state solution, we use 
pscudo-tinv; steps with a fixed Courant 
number. For time dcpicndcnt calcula- 
tions, for reasons of both stability and 
efficiency, it is best to take several 
smaller time steps on the fine grids for 
every one coarse grid step. We have 
experimented with taking several steps on 
a large fine grid for every coarse grid 
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Figure 2.3 Two fine grids generated 
around one group of flagged points. 

step. The optimal strategy for the 
number of iterations to be done on each 
grid at a time is still an open question. 
This will be an especially important con- 
sideration in large computations where 
secondary storage is used. 

We briefly review the finite volume 
Runge Kutta time stepping procedure 
which is the integrator for this adaptive 
algorithm. For details, see [12] and the 
discussion of the FL052 program. The 
full Euler equations in 2D are written in 
integral form, 

■—ffwdxdy + f^^dy-gdx = 0, 
where 


p 

PM 

/ = 

P« 

PM^+P 

g = 

pv 1 
1 

pv 

pMV 


IP^J 

1 

J 


pvH J 


The equations arc approximated in a com- 
putational domain where the variables are 
ocU-ocntcred. The flux is evaluated at the 
boundary of a given cell using the aver- 
age of the values in the adjacent cells. 
l%ds spatial discretization procedure leads 
to a system of ordinary differential equa- 
tions. These have the form 

^{hw) 4- Qw-" Dw - 0, 
where Qw is the aproximation to the 


Euler terms, Dw is an added dissi{»tive 
term, and where b is the cell area. The 
dissipation is introduced by a combination 
of second and fourth orwr differences, 
which are switched on by pressure gra- 
dients. The same dissipation formulas are 
used in the integration step on each grid, 
except at the boundaries of the fine grids, 
where the fourth order stencil is too 
large. In this case we use only the second 
order dissipation. 


The ODE'S arc Integrated using a 
modified four stage Runge Kutta scheme 
in which the dissipative terms arc only 
evaluated once. At each time step the 
solution is updated by the following 
sequence; 


«,(!) 

vv(2) 


W'~' 


= H»(0) 


/(\\ 

W'~' 


w 


(4) = ^(0) _ 


At fo\ ^ mv 


where is the value at the beginning 
of the time step, and is tte fin^ 
updated value. The time step limit for 
this scheme is almost the same as that of 
a standard fourth order Runge Kutta 
scheme. 


The farficld boundary values are 
partially specified from frecstream values, 
and partially extrapolated from the 
Ricmann invariants, depending on 
whether the flow is supersonic or sub- 
sonic, and whether the boundary is an 
inflcw or outflow boundary. At the 
brdy, only the pressure is needed to 
aJvphoe the variables in the cells nearest 
tLi ^ ill. The pressure is computed by an 
extt .folation formula based on the nor- 
mal momentum equations. The slight 
F odification of this required at the 
cnarsc/finc interface is dc^bed in sec- 
tion 3. 


3. Fine Grid Boundary Conditions 

In this section we discuss the differ- 
ence equations used at the interface 
between a coarse and fine grid. The pro- 
oedure which we have developed is 
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designed for use with the finite volume 
difference schen>c. More general pro- 
cedures are described in (13J. 



Flmra 34 Fine/ coarse arid inter- 
face. 

We illustrate the procedure using a 
refinement ratio of 2 between grids. In 
figure 3.1, the coarse grid variables arc 
marked with o, the fine grid variables 
with an x. Notice that there arc coarse 
points "underneath’* the fine grid. The 
solution is computed over the entire 
coarse grid including these points when a 
coarse grid integration step is taken. 
However, the coarse grid [mints under- 
neath the fine grid arc updated after each 
coarse step by replacing the solution there 
with the volume weighted average of the 
solution at the four nearest fine points. 
Tbe final output uses the values from the 
finest grid covering each region. 

To take a step on the fine grid, the 
flux across the line 4=0 into the fine grid 
must be calculated. It is algorithmi^y 
advantageous to introduce an "outside" 
column of solution values (denoted by the 
plus signs in figure 3.1). The flux can 
then be calculated in tbe same way for 
the first cell and the interior cells of the 
fine grid. This also mimics the coarse 
grid seti;p, where an extra column is kept 
on each side of the periodic boundary. 
The easiest way to determine this column 


of fine grid values is by interpolation 
from the coarse grid. Point p would be 
set by linear interpolation from coarse 
grid points To 

maintain accuracy, the values at V/+ij 
and V/ 4 . 1 /+! on the coarse grid would be 
replaced oy the volunw weighted average 
of the four neighboring fine grid points 
after every coarse grid integration step. 
Jb this way, each grid can still be 
integrated independently, in a manner 
that still vectorizes, and only a small 
amount of "fixup" work along the boun- 
daries of the fine grids need be done. 


Unfortunately, there is no reason for 
this procedure to be conservative, which 
in this case means that the sum of the 
fluxes into the coarse cells at the interface 
(computed for example using the value 

V/+W 

2 ) equal to that computed 

on the fine grid into the fine cells on the 
right of the interface. It is important to 
maintain conservation in order to guaran- 
tee the correct shock location in transonic 
flow fields. This is especially relevant 
since there will often be a fuic grid in the 
region of a shock, and so the interface 
between the fine and coarse grids will be 
near the shock. 


An alternative interface procedure 
which is conservative is to ailculatc the 
flux on ^hc coarse grid, and divide it in 
half for the adjacent two fine cells. This 
would bypass setting the outside vari- 
ables, and calculate a boundary flux for 
the fine grid directly. Unfortunately, this 
is unstable, as can be seen from a linear 
analysis of the interface. The treatment 
in tbds ease yields the linear relationship 

Vo,l+Vi,l + Vo,2+Vi 

where u approximates the solution on the 
coarse grid, and v approximates the solu- 
tion on the fine grid. It turns out that 
any boundary scheme which couples the 
Unc points across the interface can give 
rise to an oscillatory wave emanating 
from the interface into the fine grid, and 
supported by the central differenced 
(Unprized finite volume) scheme. 
Another alternative, that of calculating 
the flux directly from the one coarse point 
and adjacent two fine points, is stable, 



but of lower order accuracy. In effect, it 
treats the solution in the column of plus 
signs as piecewise constant in each coarse 
cell, instead of linear. 

The procedure we have chosen is a 
variation of interpolation. The cells with 
plus signs arc obtained from interpolation 
from ±e coarse grid, and die fine grids 
arc then advanced, The coarse grid 
fluxes arc determined as usual during the 
regular integration step, but then the 
value at c^ch coarse grid point nearest the 
interface is "fixed" so that the flux at the 
interface equals the sum of the two fine 
cell fluxes. In this way, conservation is 
enforced, second order accuracy is mainv' 
tained, and the only extra computational 
work is done along the boundary. (TTiis 
method thus avoids having to chc^ every 
coarse point to determine whether it is 
located at an interface, which would be 
unacceptable overhead). In experiments 
with these interface equations, when tk‘ 
interface is forced to be right at the locsi* 
tion of a shock, no loss of accuracy is 
observed in the solution. 

A special treatment is required when 
the interface coincides with the body. In 
order to interpolate for a fine point value 
at the body, a coarse grid value is 
needed at the body as well. Since only the 
pressure is computed at the body, values 
arc needed for the density, x and y 
momentum, and energy. These are com- 
puted by setting the momentum normal to 
the body to 0, setting the tangential 
momentum to be identical to that one cell 
over, and setting the energy to its steady 
state constant vdue. In practice, there is 
little difference between this procedure 
and simple linear extrapolation of the 
missing coarse grid values from the inte- 
rior. 

4. Numerical Results 

We present results comparing the 
solution computed on a coarse grid, on a 
coarse grid with patched fine grids, and 
on a ui^ormly refined grid. In all eases, 
by refining a fraction of the grid, the 
accuracy of the solution on a uniformly 
fine grid is recovered at less than half the 
cost. 


The first test ease is for non-lifting 
subsonic flow over a NACA 0012 airfoil. 
The Mach number is .500 with xcro 
degrees angle of attack. Figure 4.1 
shows a grid ^cf 'tkm 32 by 8, with the cal- 
culated pressure coefficient shown in Fig- 
ure 4.2. The drag coefficient is .0049. 
Figure 4.3 shows a grid of size 64 by 16, 
'With the pressure coefficient in Figure 
4.4, The drag coefficient is .0011. Fig- 
ure 4.6 is the solute on a grid of si^ 
128 by 32 grid, shown in Figure 4.5. The 
drag coefficient here Is ,0W2. TJic drag 
coefficient is converging like hr to its 
expected value of zero. Figure 4.8 shows 
a refined grid solution based on a 32 8 

underlying coarse grid, with refined grid 
patchcf^ as shown in Figure 4,7. The drag 
cocfficilent in this ease is .0009. Figure 

4.10 shows a refined grid solution based 

on a 64 by 16 underlying coarse grid. 
The refined grid for this case is shown in 
Figure 4.9, drag coefficient is 

reduced to .UUOl. In boQi cases, file 
accuracy of the solution on the uniformly 
next finer level grid is recovered by using 
small grid patches at the leading and trail- 
ing edges of the airfoil. 

The second test case is transonic 
flow containing a shock wave. Figure 

4.11 shows the pressure coefficknt for a 
NACA 0012 airfoil at Mach .8 with zero 
degrees angle of attack. The mesh used 
for this computation is 64 by 16 cells. 
When the grid is refined (using an error 
tolerance of .005), as shown in Figure 
4.12, the solution obtained is almost 
identical to the solution computed on a 
128 by 32 mesh (compare Figures 4.13 
and 4.14). In the coarse grid run, the 
entropy behind the shock was computed 
to be .0072, in the multiple grid run it 
was .0052, and in the fine grid run the 
entropy was .0054. In this mesh refined 
solution, 21 % of the coarse grid was 
refined by a factor of 2 in both coordinate 
directions. The cost of integrating the 
mesh refined run was thus roughly half 
the cost of the 128 by 32 grid run. If the 
error toleraiioe for m^ refinement is less 
stringent (.025), so that only the leading 
and trailing edges are refined (Figure 
4.15), the solution is only slightly worse 
(Figure 4.16). The entropy pr^ction 
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across the shock is .0046 in this ease. In 
this run only 10 % of the coarse grid is 
refined, and so the ovcmil cost of the 
solution is roughly 35 % of the flne grid 
cost. 
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